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HORIZON-BASED RESIDUAL DEPTH MIGRATION 

VELOCITY ANALYSIS 

Cross-reference to related applications 

Not applicable. 

Statement regarding federally sponsored research or development 

Not applicable. 

Background of the Invention 

Field of the Invention 

[0001] The invention relates generally to the field of seismic data processing. More 

specifically, the invention relates to methods for migrating seismic data to correct arrival 
times and apparent depths of reflective events for changes in the contour of subsurface 
reflective horizons and for changes in velocity of the formations through which seismic 
energy propagates. 

Background Art 

[0002] Seismic surveying is used to determine structures of, to determine compositions 

of, and to determine fluid content of subsurface earth formations, among other uses. A 
particular application for seismic surveying is to infer the presence of useful materials, 
such as petroleum, in the subsurface earth formations. Generally, seismic surveying 
includes deploying an array of seismic sensors at or near the earth's surface at selected 
geographic positions, and deploying one or more seismic energy sources at selected 
locations, also at or near the earth's surface. The one or more seismic energy sources are 
actuated and seismic energy emanates from the source(s), traveling generally 
downwardly through the earth's subsurface until it reaches one or more acoustic 
impedance boundaries in the earth. Seismic energy is reflected from the one or more 



1 



PGS-03-07US 



impedance boundaries, where it then travels upwardly until being detected by one or 
more of the seismic sensors. Structure and composition of the earth's subsurface is 
inferred from the travel time of the reflected seismic energy, from the geographic position 
of the source to each of the sensors and from the amplitude and phase of the various 
frequency components of the reflected seismic energy with respect to the energy 
emanating from the seismic source. 

[0003] Structures of the earth's subsurface are inferred from the travel time of the 

seismic energy from the source to the acoustic impedance boundaries and back to the 
seismic sensors at the surface. In order to infer depth of and the structures of subsurface 
earth formations from reflection seismic travel times measured at the earth's surface, it is 
necessary to determine the acoustic velocity of the various formations through which the 
seismic energy passes. Velocities of the earth formations can vary both with respect to 
depth in the earth (vertically), and with respect to geographic position (laterally). 
Seismic data, however, are recorded only with respect to time. Methods known in the art 
for estimating velocities of the earth formations both vertically and laterally rely on 
inferences about the travel path geometry of the seismic energy as it travels from the 
seismic source to the various seismic receivers deployed at or near the earth's surface. 

[0004] Migration is a process performed on seismic data in which depth estimates to one 

or more reflective horizons (acoustic impedance boundaries) in the earth are made from 
the "two-way" travel time of seismic energy from the seismic energy source /to the 
reflective horizons and back to the seismic receivers. The depth estimates of the 
reflective horizons are computed and are displayed with respect to geographic position. 
Depth estimates based on two-way travel time must be corrected for energy travel path 
differences between the various seismic energy source and receiver geographic positions 
that are used during data acquisition. In order to correct the depth estimates for the 
various source and receiver positions, it is necessary to accurately estimate the velocity of 
seismic energy in the earth from the earth's surface (or the ocean bottom in marine 
seismic data) to each subsurface reflective horizon. Methods are known in the art for 
estimating velocity from two-way travel time from the seismic source to the reflective 
horizons and back to the seismic receivers. One such method uses two-way travel times 

2 



PGS-03-07US 



for source and receiver arrangements which have a "common mid point" along the 
seismic energy travel path. Acoustic velocities of the earth formations from the earth's 
surface to a particular subsurface reflector can be estimated using the familiar Dix 
equation, for example. Other methods for estimating velocity are known in the art. 

[0005] According to wave propagation theory well known to those skilled in the art, a 

spherical seismic energy wave propagating from a "point" source (a source modeled for 
calculation purposes as having essentially zero volume or spatial extent) can be 
decomposed into a series of plane waves. See, for example, Stoffa, P.L., Buhl, P., 
Diebold, J.B., and Wenzel, F., Direct mapping of seismic data to the domain of intercept 
time and ray parameter — A plane-wave decomposition: Geophysics, 46, 410-421 (1981). 
The Stoffa et al. article describes a method for decomposing seismic reflection data into 
plane waves by slant stacking, i.e., transforming seismic data into the plane wave (x-p) 
domain (x = intercept time and p = ray parameter), and further documents some 
advantages of processing seismic data in the plane wave (r-p) domain, including, for 
example, linear noise attenuation or normal move-out in a horizontally stratified medium 
without approximation. Using stacking velocity analysis performed in the offset-time (x- 
t) domain, by contrast, provides RMS (root mean square) velocities. The RMS velocities 
can then be converted into interval velocities when required. Advantageously, velocity 
analysis in the plane-wave domain results in the estimation of interval velocities directly. 
Having accurate estimates of interval velocities is important for performing migration. 

[0006] Some of the research in prestack migration velocity analysis began in the early 

1990's. See, for example, Al-Yahya, K., Velocity analysis by iterative profile migration, 
Geophysics, vol. 54, pp. 718-729 (1989). Various analytic functions have been derived 
to express the relationship between the true velocity (or the ratio of the migration velocity 
and the true velocity) and the offset in a common image gather (CIG) in the depth-offset 
domain after migration. The foregoing analytic functions make use of the assumptions of 
a small dip (rate of change of depth with respect to lateral displacement), small offset, 
and/or constant velocity in the various layers of the earth's subsurface. Residual moveout 
analysis has also been used to extend the application of such analytic functions to media 
having lateral velocity variation. See, for example, Meng, Z, Bleistein, N, and Wyatt, 
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K.D, 3-D Analytical migration velocity analysis I: Two-step velocity estimation by 
reflector-normal update, 69 th Annual International Meeting, Society of Exploration 
Geophysicists, Expanded Abstracts (1999). 

[0007] Most of the migration methods known in the art are implemented in the depth- 

offset domain (z-x\ and a top-down "layer stripping" migration method is then used to 
derive the interval velocities. It is known in the art to use the depth-offset domain 
because this is the domain in which most prestack depth migration is performed, and the 
domain in which migrated CIG's are available for analysis. However, it is also known in 
the art to perform prestack depth migration in the plane- wave (z-p) domain. See, for 
example, Akbar, F.E., Sen, M.K., and Stoffa, P.L, Prestack plane-wave Kirchhoff 
migration in laterally varying media, Geophysics, 61, 1068-1079 (1996). See also, 
Tanis, M. C, Prestack Split-step Fourier Depth Migration Algorithms and Parallel 
Implementation on Cray T3E, Ph.D. Dissertation, The University of Texas at Austin 
(1998). After migration in the plane wave domain, seismic data are displayed or 
presented in the depth-plane wave (z-p) domain. Prestack depth migration using slant 
stack (z-p) data and a substantially correct interval velocity-depth model generate events 
in a common image gather (CIG) in the depth-plane wave (z-p) domain which are 
substantially horizontally aligned, because a CIG represents an image of the same 
subsurface position obtained along different seismic travel path angles. See, for example, 
Whitmore, N. D. and Garing, J. D., Interval velocity estimation using iterative prestack 
depth migration in the constant angle domain, The Leading Edge, vol. 12, no. 7, pp. 757- 
762 (1993). 

[0008] Use of an erroneous velocity-depth model in migration, however, can cause 

misalignment of reflective events in a CIG, meaning that the reflective events displayed 
on the CIG exhibit a residual "moveout." By analyzing the residual moveout (a change 
in apparent depth with respect to ray parameter) in the CIG, it is possible to derive depth 
and velocity corrections, thus obtaining an updated velocity-depth model. For example, 
if the velocity used in the migration process is lower than the true velocity, the event 
appears to curve upwardly in the depth-plane wave (z-p) domain after prestack depth 
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migration. If the velocity used in the migration process is higher than the true velocity, 
then the events in the CIG appear to curve downwardly. 

[0009] For some time, a method known as the "vertical velocity update method" has been 

used to generate a velocity-depth model for prestack depth migration. A typical data 
processing procedure used in such methods is known as the "Deregowski loop." See 
Deregowski, S. M., Common-offset migrations and velocity analysis, First Break, vol. 8, 
no. 6, pp. 224-234 (1990). Residual velocity analysis can be applied at all depths based 
on the constant velocity assumption. See, Al-Yahya, K. (1989), Velocity analysis by 
iterative profile migration, Geophysics, vol. 54, pp. 718-729. Then the constant 
velocities are converted to interval velocities for a subsequent iteration. If it is desired to 
obtain the interval velocities from migrated seismic data directly, it is necessary to 
perform both prestack depth migration and the velocity analysis in a top-down "layer- 
stripping" procedure. 

[0010] More recently, a method has been developed to update interval velocities using 

residual analysis in the depth-plane wave domain. See, Jiao, J., Stoffa, P., Sen, M., and 
Seifoullaev, R., Residual migration velocity analysis in the plane-wave domain, 
Geophysics, vol. 67, pp. 1258-1269 (2002). See also, Jiao, J., Residual Migration 
Velocity Analysis in The Plane Wave Domain: Theory and Applications, Ph.D. 
Dissertation, The University of Texas at Austin (2001). The method disclosed in the 
foregoing reference eliminates the need to perform "layer-stripping" prestack depth 
migration in order to obtain interval velocities. However, the method disclosed in the 
Jiao et al. article requires that prestack migration be performed in the depth-plane wave 
domain (z-p), which limits the application of that method. It is desirable to have a 
method for performing prestack migration in the depth-offset (z-p) domain which 
includes updating of interval velocities. 

Summary of the Invention 

[0011] One aspect of the invention is a method for processing seismic data. The method 

includes prestack depth migrating the seismic data using an initial velocity-depth model. 
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At least one horizon in the migrated seismic data is selected. Residual migration velocity 
analysis in the depth-offset domain is performed with respect to the at least one selected 
horizon, and the velocity-depth model is updated based on the residual migration velocity 
analysis. In some embodiments, a deeper horizon is selected, and residual migration 
velocity analysis is again performed on the deeper horizon. Based on the residual 
migration velocity analysis on the deeper horizon, the velocity-depth model is again 
updated. 

[0012] Another aspect of the invention is a computer program stored in a computer 

readable medium. The program contains logic operable to cause a programmable 
computer to perform the following steps. Seismic data are prestack depth migrated using 
an initial velocity-depth model. At least one horizon in the migrated seismic data is 
selected. Residual migration velocity analysis in the depth-offset domain is performed 
with respect to the at least one selected horizon, and the velocity-depth model is updated 
based on the residual migration velocity analysis. In some embodiments, a deeper 
horizon is selected, and residual migration velocity analysis is again performed on the 
deeper horizon. Based on the residual migration velocity analysis on the deeper horizon, 
the velocity-depth model is again updated. 

[0013] Other aspects and advantages of the invention will be apparent from the 

description and claims which follow. 

Brief Description of the Drawings 

[0014] Figure 1 shows a flow chart of a prior art residual migration velocity analysis 

method. 

[0015] Figure 2 shows a flow chart of an embodiment of a method of the invention. 

[0016] Figures 3a and 3b show, respectively, a depth migrated seismic section using a 

velocity depth model from a conventional velocity analysis method, and a depth migrated 
section using a velocity depth model derived from residual migration velocity analysis 
according to the invention. 
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[0017] Figure 4 shows a programmable computer used to read a program stored in a 

computer readable medium. 

Detailed Description 

[0018] For purposes of understanding the invention, it is helpful to review briefly the 

process of residual migration velocity analysis in the depth-plane wave domain. A trial 
slowness (the inverse of velocity), represented by a/ 1 , of the /* layer of a substantially 
horizontal, plane-layered medium, having a total number N of such layers, is used for 
prestack depth migration. Any difference between the true slowness uj of the /* layer 
and the trial slowness uf will result in a misfit between the true depth z and the 
migration-calculated depth z m . The misfit is referred to as the residual difference, z reS9 or 
"residual moveout" and can be determined with respect to the ray parameter p by the 
following expression: 

[ooi9] z r jp)=±Azj(py 

[0020] where uj and Azj represent, respectively, the true slowness and the vertical 

thickness of the/ layer, and Az™ represents a trial thickness of the/ layer. See, Jiao, J., 
Stoffa, P., Sen, M., and Seifoullaev, R., Residual migration velocity analysis in the plane- 
wave domain, Geophysics, 67, 1258-1269 (2002) for a derivation of the foregoing 
residual moveout equation. In another implementation, the residual difference z res (p) can 
be calculated by the expression below: 

[00211 ».W - j>7(|.-0).i£.(l- - "']" I (la) 

y-i w, [ \\Mj ) -P J J 

[0022] Equation (la) is more efficient in some practical applications than equation (1) 

because equation (la) contemplates calculation of the depth residual using only the layer 
thickness corresponding to a ray parameter value equal to zero. 
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[0023] Equations (1) and (la) each describe the depth residuals in depth-migrated results 

as a function of the ray parameter /?, the true slowness u } and the trial slowness wj . The 

true depth, z, of a particular subsurface earth layer can be calculated from the migration 
calculated depth and the residual by the following equation: 

[0024] z = z m (p) + z res (p) (2) 

[0025] A practical implementation of the process described with reference to equations 

(1), (la) and (2) used to perform depth migration in the depth-plane wave domain is 
shown as a flow chart in Figure 1. At 10, seismic data are recorded with respect to the 
time of actuation of a seismic energy source. The data recordings are each identified with 
respect to the individual receiver which made the particular recording. Each recording 
thus has associated therewith a distance between the source and the individual receiver 
(called "offset") at the time of recording. Thus, the recordings may be mathematically 
represented in the time-offset (x-t) domain. At 12, the recordings are transformed from 
the time-offset domain into the plane wave (x-p) domain. An initial velocity model, at 
16, is used to perform prestack depth migration, at 14, in the plane wave domain. At 18, 
common image gathers ("CIGs") of the migrated seismic data are generated. At 20, one 
or more of the CIGs are selected for processing, and at 22, residual migration velocity 
analysis is performed as explained above with respect to equations (1) or (la), and (2) on 
the selected CIGs for a first selected reflective event or "horizon." The migration 
residual velocity analysis includes calculating the depth residual with respect to the ray 
parameter using either equation (1) or (la). A small perturbation within a selected range 
is added to the migration slowness, and the depth residual with respect to ray parameter is 
calculated again. The perturbation is changed and the depth residual calculation is 
repeated until the entire selected range of perturbations has been applied to the migration 
velocity. The perturbation which results in the "flattest" apparent depth of the first 
reflective horizon is then selected. The selected perturbation is then added to the trial 
slowness to generate an estimated true slowness for the first reflective horizon. 

[0026] At 24, the initial velocity model and a layer thickness model are adjusted with 

respect to the residual migration velocity analysis performed at 22. At 26, if additional 
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reflective horizons deeper than the first reflective horizon are to be processed, the 
selected CIGs are then subjected to a residual velocity analysis, at 22, with respect to the 
next selected event or horizon in the each of the CIGs. As with the first selected event or 
horizon, in the next, and in any successive horizons, a range of perturbations is applied to 
the slowness of the particular layer, and the perturbation which when added to the initial 
slowness for the layer provides the flattest relationship of migration depth with respect to 
ray parameter is added to the initial slowness to estimate the true slowness for the 
particular layer. 

[0027] At 28, if additional CIGs are to be processed, the process returns to CIG selection, 

at 20. If all selected CIGs have been processed, then the velocity-depth model is updated, 
at 16, and prestack depth migration can be performed again, as shown at 14, using the 
updated velocity-depth model. The process can be repeated until all the selected CIGs 
are processed. At 30, the updated velocity-depth models can be interpolated into an 
entire line of seismic recordings. 

[0028] The foregoing process can be extended from the depth-plane wave domain to the 

depth-offset domain. The theoretical basis for such extension can be understood as 
follows. For a seismic energy ray traveling from a seismic energy source to a seismic 
receiver through a substantially horizontally stratified medium, the following relationship 
exists between the offset x and the ray parameter p: 



represents the distance between the seismic source and the particular receiver (the offset) 
and p represents the ray parameter. The ray parameter p can be determined with respect 
to the offset x by virtually inverting equation (3): 



[0029] 




(3) 



[0030] 



in which V(z) represents the velocity of the medium with respect to depth, x 



[0031] 



p(x) = inv 



{*(/>)} 



(4) 
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[0032] Substituting equation (4) into the depth residual equation (1), the depth residual 

for a 1 -dimensional medium in the depth-offset domain can be determined from the 
slowness and the offset x by the expression: 



[0033] Z res (x,p(x))=fl 

7=1 



&Z"(xi J J — l 

' 1 iu)-p\x)r 



(5) 



[0034] in which u Jt uf and p represent the same parameters as in equation (1). If depth 

residual equation (la) is used for the residual depth analysis instead of equation (1), 
equation 5 may be substituted by the following: 



[0035] Z res (x) = XAz?(?c = 0).^.{l- 



(5a) 



[0036] After determining the residual depth correction, a corrected depth in the depth- 

offset domain can be determined for each subsurface layer using the expression: 

[0037] z(x) = z m (x) + z res (x) (6) 

[0038] The residual depth correction determined using equation (6) can then be used for 

migration velocity analysis. 

[0039] In deriving the depth residual equations above, it was assumed that both the 

migration slowness uf 1 and the true slowness uj are known for each layer. As a practical 
matter, however, the true slowness uj is the parameter to be determined by migration 
velocity analysis, and the migration slowness uj m is an approximation of the true slowness 
used when performing the migration. Thus, as in the case for residual depth analysis in 
the plane wave domain, residual depth analysis in the offset domain can be performed by 
applying a perturbation from within a range of perturbations to the migration slowness 
u/ 1 to generate a new slowness value uj new = up + up\ In the foregoing expression, u es 
represents the slowness perturbation and may be called the "residual slowness." The 
closer the new slowness w/ evv is to the true slowness, the "flatter" reflective events will 
appear on a common image gather (CIG) with respect to offset, after applying the 
residual depth corrections. In one implementation, an initial value of slowness 
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perturbation is applied to the initial value of migration slowness, and the flatness in the 
CIG is determined with respect to offset along a first selected event or horizon. The 
slowness perturbation value is changed, and a new value of migration slowness is then 
determined. CIG event flatness is again determined. Changing the slowness perturbation 
value, and calculating the CIG event flatness is repeated until the entire selected range of 
slowness perturbation values is covered. The value of the slowness perturbation which 
results in the greatest degree of flatness is added to the migration slowness value to 
provide a closest approximation of the true slowness for the particular layer being 
evaluated. The velocities can be determined as the inverse of the true slowness values 
thus determined. The foregoing process can then be repeated for any additional 
subsurface horizons desired to be analyzed in each CIG. 

[0040] To obtain interval velocities using the above process, a top-down "layer stripping" 

technique is used. After performing the foregoing residual velocity analysis along a first 
selected horizon (layer) based on the initial velocity model used for the depth migration, 
velocities V" for the first layer are updated to V 1 ™. Depths 2T (or thickness) for the first 
layer are updated to Z" ew by following equation: 

Trnew 

[0041] Z =Z — (7) 

[0042] Residual migration velocity analysis can be performed as explained above for the 

next selected horizon downward in depth, based on the updated velocity-depth model. 
The process can then be repeated for each selected horizon until all selected horizons are 
processed. After completing the foregoing process for all horizons at any CIG, a final 
velocity-depth model is obtained that is typically closer to the true velocity field than the 
previous model iteration or the initial model used in the depth migration. The final 
velocity-depth model can be used for a subsequent iteration of prestack depth migration, 
or can be used only to perform a residual depth or velocity correction on any one or more 
CIGs. 

[0043] Determining the "flatness" in a CIG may be performed by determining semblance 

between traces in the CIG. When semblance reaches a maximum for a selected value of 
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linear depth shift with respect to offsets, the CIG is determined to be at maximum 
flatness. In some embodiments, a selected threshold value for semblance (flatness) may 
be used instead of a maximum value. 

[0044] Figure 2 shows a flow chart of one implementation of a velocity analysis method 

of the invention used in seismic data processing. At 10, seismic data recorded with 
respect to time and receiver location (as explained above with respect to Figure 1) are 
selected for processing. At 32, prestack depth migration can be performed, for example, 
using the Kirchhoff integral method. The prestack depth migration performed at 32 
includes an initial velocity-depth model, at 42. At 33, common image gathers are 
generated in the depth offset domain using the velocity model from the prestack depth 
migration. At 34, a first reflective event, or horizon, is selected for processing. At 35, 
residual velocity analysis is performed in the depth-offset domain. The residual velocity 
analysis, as previously explained, can be performed by applying a perturbation within a 
selected range to the velocity from the initial migration (at 32), and determining the 
flatness of the selected horizon in the CIG. The velocity perturbation which when added 
to the initial value of velocity results in the greatest degree of flatness is added to the 
initial value of velocity to determine the true velocity for the layer being processed. At 
37, the layer velocity, layer depth and layer thickness are updated based on the residual 
velocity analysis performed previously at 35. At 39, the residual velocity analysis is 
performed along successively deeper horizons, as explained above with respect to 
equations (5) through (7), until the entire velocity-depth model is updated, at 36. At 38, 
the updated velocities determined from the residual analysis may be used to perform 
another depth migration. The new depth migration will result, at 42, in an updated 
velocity-depth model. The process can be repeated until all selected horizons have been 
processed, or at 40, to perform residual corrections on common image gathers. 

[0045] Methods according to the invention extend residual normal moveout (NMO) 

velocity and depth correction from the depth-plane-wave domain to the depth-offset 
domain and provide a residual depth migration velocity analysis method using common 
image gathers in the depth offset domain. To obtain interval velocities using methods 
according to the invention, it is only necessary to perform top-down layer-stripping 
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residual corrections instead of top-down layer-stripping prestack migration. Thus 
methods according to the invention are particularly suitable for migration-velocity 
analysis, especially for large three-dimensional seismic surveys. 

[0046] In a three-dimensional field data test, a good velocity-depth model was obtained 

after only two iterations. The velocity-depth model obtained using the field test data are 
shown in Figure 3a, using prestack depth migration such as the Kirchoff integral method. 
Residual migration velocity analysis in the offset domain was then used to migrate the 
entire three-dimensional data volume again. The image was improved significantly as 
compared to the original depth migration, as can be observed in Figure 3b. Methods 
according to the invention can also be used to perform residual correction on an entire 
data set that will improve the quality of stacked common image gathers. 

[0047] In another aspect, the invention relates to computer programs stored in a computer 

readable medium. Referring to Figure 4, the foregoing process as explained with respect 
to Figure 2, can be embodied in computer-readable code stored on a computer readable 
medium, such as floppy disk 64, CD-ROM 62 or magnetic hard drive 606 forming part of 
a general purpose programmable computer. The computer, as known in the art, includes 
a central processing unit 50, a user input device such as a keyboard 54 and a user display 
52 such as a flat panel LCD display or cathode ray tube display. According to this aspect 
of the invention, the computer readable medium includes logic operable to cause the 
computer to execute steps as set forth above and explained with respect to Figure 2. 

[0048] While the invention has been described with respect to a limited number of 

embodiments, those skilled in the art, having benefit of this disclosure, will appreciate 
that other embodiments can be devised which do not depart from the scope of the 
invention as disclosed herein. Accordingly, the scope of the invention should be limited 
only by the attached claims. 
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